Ergodic Theory and Visualization I: Mesochronic Plots for 
Visualization of Ergodic Partition and Invariant Sets 
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We present a computational study of a visualization method for invariant sets based on ergodic partition 
theory, first proposed in [TJ|2]. The algorithms for computation of the time averages of observables on phase 
space are developed and used to provide an approximation of the ergodic partition of the phase space. 
We term the graphical representation of this approximation - based on time averages of observables - a 
Mesochronic PlotQ The method is useful for identifying low-dimensional projections (e.g. two-dimensional 
slices) of invariant structures in phase spaces of dimensionality bigger than two. We also introduce the 
concept of the ergodic quotient space, obtained by assigning a point to every ergodic set, and provide 
an embedding method whose graphical representation we call the Mesochronic Scatter Plot (MSP). We 
use the Chirikov standard map as a well-known and dynamically rich example in order to illustrate the 
implementation of our methods. In addition, we expose applications to other higher dimensional maps 
such as the Froeschle map for which we utilize our methods to analyze merging of resonances and, the 
three-dimensional Extended standard map for which we study the conjecture on its ergodicity [SJ. We 
extend the study in our next paper [3] by investigating the visualization of periodic sets using harmonic 
time averages. Both of these methods are related to eigenspace structure of the Koopman operator. 



"From Greek: meso - mean, chrono - time 

Dynamical equations describing behavior of systems 
of scientific interest are often impossible to solve 
analytically, and one must resort to various com- 
putational methods approximating the actual solu- 
tion. We develop computational methods for study- 
ing invariant sets of dynamical systems with arbi- 
trary dimensionality using methods of ergodic par- 
tition theory. Our method consists in computing 
time averages of chosen functions under the phase 
space dynamics for a given trajectory. This allows 
us to obtain graphical visualization of invariant sets, 
including easily obtainable intersections of such sets 
with lower-dimensional surfaces or manifolds. The 
numerical efficiency of the method - interestingly - 
improves when the dynamics is more complex. To 
represent the obtained information in a condensed 
form, we introduce scatter plots that represent the 
global topology of ergodic sets. We illustrate the 
utility and extent of our methods by investigating 
the global dynamical properties of several known 
measure-preserving dynamical systems. In the con- 
text of dissipative systems, the method visualizes 
basins of attraction. 



1 Introduction 

The increasing range and dimensionality of phenomena mod- 
eled as dynamical systems creates the demand for new and 
better computational approaches to comprehend high-dimensional, 
complex dynamics. Most of the dynamical systems prob- 
lems of current scientific interest cannot be treated by purely 
analytical techniques, constraining one to chose among di- 
verse computational methods [5]. The choice of the numer- 
ical method is primarily determined by the nature of the 
problem, but also by the type of investigation that is to be 
conducted - the solution of the entire problem is often not 
required, and the attention is focused on the relevant results 
only. 

The choice of the computational approach can be the 
decisive factor for the overall efficiency of investigation and 
the precision of the results. A variety of methods are used 
for studying dynamical systems - the direct numerical inte- 
gration of equations of motion still remains the common ap- 
proach, both for discrete-time (maps) and continuous-time 
(ODEs) systems. In this context, the visualization methods 
play a major role: it is often of interest to graphically visu- 
alize various aspects of motion, for instance the phase space 
structure of a given dynamics, without necessarily solving 
the entire system of equations. However, there is lack of 
such methods that easily extend to high-dimensional, com- 
plex nature of motion exhibited by many dynamical sys- 
tems of interest. It would be useful, in particular, to have a 
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method that enables visualization of low-dimensional (e.g. 
2D) slices through higher dimensional invariant structures. 
We develop such a method, suggested first in [TJE]. 

Visualization methods generally work as algorithms for 
dividing the dynamics' phase space into subsets according 
to a prescribed property of interest. Exit time plots [5] are 
computed by fixing a bounded subset of the phase space and 
measuring the time needed for its representative points to 
exit the subset, which are colored according to their exit- 
times. The phase space is then sliced in the equal-exit-time 
regions, giving conclusions regarding the system's transport 
properties. The set-oriented computational algorithms [7] 
allow visualization of the invariant sets and attractors by 
the appropriate phase space subdivision. After covering a 
tentative invariant region with boxes of certain (small) size, 
a sequence of box-size reductions is applied, and the opti- 
mal approximation of the invariant geometry is obtained at 
the limit. Similarly, by studying the return time dynam- 
ics [8] one constructs a phase space decomposition into al- 
most invariant sets by measuring and comparing the return 
times for each subdivision element. Methods of computing 
(un)stable manifolds based on analysis of geodesic level sets 
were proposed in the context of vector fields [S], having a 
particular importance for control problems. Invariant man- 
ifolds can also be visualized using distribution of finite-time 
Lyapunov exponents in phase space [10 , or by integration 
of fat trajectories which extend the concept of standard nu- 
merical trajectory integration [llj . 

On the other hand, ergodic theory [12] as a study of sta- 
tistical aspects of motion is extensively used in the context of 
chaotic dynamical systems [T3] • Its applications gave results 
that range from optimal mixing |14j . traffic jams |15j and 
quantum many-body problems [IB], to mathematical study 
of maximizing measures in discrete dynamics [17] . Finally, 
study of properties of time averages and harmonic averages 
of dynamical systems gave rise to new methods for visualiza- 
tion of invariant and periodic structures, both theoretically 
[H [TH] and experimentally [TH] • 

This paper presents a detailed computational study of 
the visualization method based on ergodic theory concepts 
proposed by Mezic in 1 . The method consists of com- 
putation of joint level sets of time averages of a basis of 
functions. This leads to an approximation to the ergodic 
partition of the phase space. We name the resulting plot 
of joint level sets the Mesochronic Plot (since the method 
is based on time averaging). This method enables study of 
high-dimensional dynamical systems by visualizing the in- 
variant set decomposition of appropriately chosen lower di- 
mensional planes/surfaces that slice/intersect the examined 
phase space region. Moreover, the method can be similarly 
employed to visualize attractor basins for systems of arbi- 
trary dimensionality. To illustrate the implementation of 
our method we start with the Chirikov standard map [5U] as 
a well-known example of a chaotic dynamical systems that 
possesses a rich variety of dynamical behaviors. We stress 
that the employment of Mesochronic Plots for 2D standard 
map is done for illustrative purposes only, while the real 
application of the method lies in systems with higher di- 
mensionality such as the extended standard map and the 



Froeschle map that we also study in this paper. 

We also introduce here the concept of the ergodic quo- 
tient space, where each ergodic set is associated with a point 
in the quotient space. We show how to approximately em- 
bed the resulting set into Euclidean space using time aver- 
ages of observables. This leads to a novel type of a plot, 
called the Mesochronic Scatter Plot (MSP). We relate prop- 
erties of MSP to dynamics in phase space. 

As opposed to the return time statistics (yielding al- 
most invariant sets) [5] and exit time [5J approaches which 
are dealing with the properties of transport and its speed, 
our method focuses on geometrical phase space properties. 
We are complementing these methods by constructing an 
algorithm that gives decomposition into ergodic sets (at the 
limit). In contrast to techniques involving finite-time Lya- 
punov exponents (similarly suitable for systems of arbitrary 
dimensionality) [10] , our method does not seek to iden- 
tify finite-time structure of normally hyperbolic invariant 
manifolds - which over time are typically densely embedded 
within the invariant structures that we define. 

In our forthcoming paper [3] which relies on theoreti- 
cal framework of Mezic and Banaszuk |18| , we employ har- 
monic time averages and develop an algorithm for visualiza- 
tion of periodic sets and resonances in the phase space. We 
show that our technique is related to known frequency map 
techniques which focus on the frequency spectrum of cho- 
sen trajectories developed by Laskar and collaborators [2"T] . 
and relate such techniques to spectral properties (in partic- 
ular, eigenfunctions) of the associated Koopman operator 
[18]. However, while Laskar et al. plot the frequencies of 
motion, we plot the phases - thus unveiling the periodic (as 
opposed to invariant) partitions in the phase space. In addi- 
tion, Laskar et al. methods are valid for near-integrable sys- 
tems, and not rigorously supported for fully non-integrablc 
systems. In contrast, our method have a rigorous justifica- 
tion for arbitrary measure-preserving systems. In relation 
to current work, while visualizing the periodic sets that res- 
onate with a particular chosen frequency is of great impor- 
tance, the technique of ergodic partitioning presented here 
provides a method of refining frequency partitioning of the 
entire phase space into even smaller invariant sets. Thus, 
frequency partitioning can in fact be thought as an example 
of a current, more general method. 

This paper is organized as follows: after briefly dis- 
cussing our method's theoretical background in Section [2j 
we show the construction of Mesochronic Plots for single 
functions under the standard map in Section [3] The con- 
vergence of time averages is addressed in Section |4j In Sec- 
tion [5] we present the ergodic quotient space concept and 
its Mesochronic Scatter Plot embedding with multiple time 
averages, and construct a simple algorithm for ergodic parti- 
tion approximation. In Sections [6] and [7] we expose concrete 
applications of our method to higher dimensional maps, like 
4D Froeschle map [23] and 3D extended standard map [3]. 
Conclusions are given in Section [8| 
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2 The Visualization Method 

In this Section we briefly sketch the mathematical basis 
of our visualization method and describe its implementa- 
tion algorithm. As our approach here is rather application- 
oriented, we will skip the rigorous proofs and refer the reader 
to pQ (and references therein) for more details. 

2.1 The ergodic theory background 

We consider a map T on a compact metric phase space 
A endowed with a measure /i that is preserved under T, 
evolving in discrete time t: 

x t+1 = Tx t , tez, x f e A. (l) 

Our central aim is to visualize invariant sets B C A for the 
map T defined as [23] : 

x e B => T'x e B Vi e Z, (2) 

which essentially means that each trajectory that starts in 
an invariant set B stays in B forever. There are many types 
of invariant sets that play substantial role in analysis of dy- 
namical systems. Here we are interested in ergodic invariant 
sets. 

Consider L 1 real-valued functions on A (/ : A — > R € 
L X (A) iff J A |/(x)|dx < oo), and let the time average /*(x ) 
of a function / 6 L 1 (A) corresponding to a phase space 
point xo £ A be defined as: 

t-i 

/*(x )=lim -]T/(T fe x ). (3) 

k=0 

By the Ergodic Theorem ([H]) this limit exists almost 
everywhere (a.e.) in A for every / € L}{A). We call the 
map T ergodic over some set B C A if for a.e. point xo € B 
we have that the time average for every function / € L 1 (B) 
is equal to the space average of that function / over the set 
B. More precisely, ergodicity of T when restricted to B (or 
ergodicity over B) means there exists an ergodic measure 
Hb such that the equality: 

f (*o) = ~4m I f d l*B, (4) 

^B\B) J B 

holds a.e. in B C A. This also implies that /* is constant 
a.e. in B C A and that almost every orbit starting in B 
covers B densely at the limit; the amount of time a tra- 
jectory spends in a given region of B is proportional to the 
ergodic measure of that region. Also, each /* is an invariant 
function /*(x ) = /*(T'xo) Vn. The set of real numbers 
K induces a partition of A called Q = {B a } ae ^ through a 
function / G L l (A) : 

B a = (f*)- 1 (a), Vaet 

where (J*) -1 is the set inverse of /*, and some of the B a 
may be empty sets. We have: 

/x((J B a ) = fi(A), B a n B a , = Va ^ a', 

a 



and, denoting the set of all the points for which the time 
average of / does not exist by £(/, T), 

A=(U a B a )|JS(/ ) T). 

While this implies each B a to be invariant by construction, 
we may still have some set B a to actually be a union of 
more independent invariant sets that accidentally have the 
same time average for a given /. For this purpose we refine 
the partitioning by considering products of partitions corre- 
sponding to different functions in order to obtain a partition 
in which each element is a non-decomposable, ergodic and 
invariant set. This final ergodic partition £ e is defined as the 
product of all belonging to a set S C ^(A): 

Ce = \/ 0' ( 5 ) 

fes 

where it suffices to take S to be a basis for L X (A) [18]. 
Defining the Koopman operator U associated with the map 
T by the composition operation 

C//(x) = / o T(x) 

it is easy to see that, the time average /* of a function / is 
an eigenfunction of U at eigenvalue 1 [18] . 

The ergodic partition has the desired properties of di- 
viding the phase space into a family of non-decomposable 
invariant sets: ergodic subsets are the "minimal observable" 
invariant sets. In order to visualize the ergodic sets we need 
to approximate the ergodic partition by computing the time 
averages of a finite number of functions, and observe the 
subsets where they are simultaneously constant (i.e. the 
joint level sets of these time averages). In the rest of this 
work we will be developing and employing a computational 
algorithm that uses the described idea. 

2.2 The computational algorithm and nu- 
merical details 

Here we describe the algorithm for approximation of the 
ergodic partition, limiting for simplicity the discussion to 
the case of a 2D map with the phase space A = [0, 1] x 
[0,1] CM 2 . 

step 1 Set up a grid (e.g. lattice) of initial grid-points 
(a^o, 2/o ) on the phase space A 

step 2 Pick N functions {/i, . . . Jn} from ^(A) and com- 
pute their partial time averages for % na i iterations for 
each grid-point, which serve as the approximations for 
the real time averages {/*, ■ - -/jy} 

step 3 To every initial grid-point (xq, yo) associate the cor- 
responding time average vector: 

(xo,y Q ) — > f (x ,yo), 

f(xo,Vo) = {/i 0eo,2/o), ■ ■•/jy(a;o,2/o)} e R w 

(6) 
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step 4 Observe the distribution of time average vectors f 
in M. N and group them optimally into clusters. Divide 
A into a union of subsets, with each subset being given 
by those grid-points (xoj2/o) whose time average vec- 
tors i(xo,yo) belong to the same cluster in R N . This 
family of subsets is an A-function approximation of 
the ergodic partition of A in the sense of Eq. |5]) . 

The optimal number of iterations i nria i is to be set in ac- 
cordance with the convergence properties of the time aver- 
ages; observe also that the quality/properties of visualiza- 
tion strongly depend on the way time average vector are 
clustered: these topic will be dealt with later. 

The choice of functions {/i,.-./iv} is important: lin- 
early dependent functions will not give new information as 
their time averages differ only by a multiplicative constant. 
It is therefore necessary to consider linearly independent 
functions, which we shall do by picking them from an or- 
thogonal basis on I/ 1 (j4). 

We use a grid of D x D initial grid-points and iterate the 
dynamics for tg na ] iterations, with values of D and ^ nna [ 
depending on the particular map/case under investigation. 

2.3 Standard map as the testing prototype 

We chose the Chirikov standard map [2U] for testing the 
performance of our method. Its behavior has been widely 
studied and well- understood [551 H3 • We consider it in the 
form: 

x' = x + y + e sin(27ra;) [mod 1] 
y' = y + esm(2Trx) [mod 1] 

where (x,y) e [0,1] x [0,1] = [0, l] 2 (the usual standard 
map's parameter k is here k — 2ire). It is an area-preserving 
(symplectic) map which exhibits a variety of invariant sets, 
both regular, composed of periodic or quasi-periodic orbits, 
and chaotic zones that evolve in size and structure as the 
parameter e is varied. 

3 Single-function Plots 

We set iV = 1 and consider the time averages of a single 
function under the dynamics of Eq.([7]) for a grid of initial 
conditions. As the range of final time average values varies, 
we adjust the coloring scheme for each plot by assigning blue 
to the minimum and red to the maximum value obtained. 

3.1 The Fourier orthogonal basis 

We pick the functions to time average from the Fourier basis 
in the form: 

foui(2n7rx) fou2(2m7ry), with n,m € N, 

where fou^ = sin or cos, for i = 1, 2, obtaining a grid of time 
averages corresponding to the grid of initial conditions. The 
time average grid is then colored, visualizing the invariant 
sets as the uniformly colored level sets. In Fig. [I] we show 
three examples of time average plots for various £-values 



(left), along with their phase space portraits (right). The 
pictures on the top deal with e = case, where the map 
possesses a family of invariant circles with y = const, (a full 
measure of these circles have ergodic dynamics - irrational 
rotation - on them) . Accordingly, the plot is a family of hor- 
izontal single-color lines, with color depending on the time 
average of the function over each y = const, line. Middle 
figures show the case of e = 0.09 where the phase space is 
a mix between a small chaotic zone and families of regular 
islands. The bottom two figures regard the case of e = 0.18 
characterized by the presence of a large chaotic zone, which 
is uniformly colored due to uniformity of time averages in 
that region. 

Utilizing a single-function partition some independent 
invariant sets appear colored with the same color in some 
plots - since the averaged functions are not one to one, same 
time average values can be obtained in dynamically differ- 
ent regions. This is why a single time average plot does 
not uniquely identify an ergodic set: this problem will be 
addressed later by employing more functions. Despite be- 
ing obtainable with usual techniques, these results nicely 
illustrate the implementation of our method. 

3.2 Multi-scale methods 

Selecting different functions to time average can reveal dy- 
namics on different scales in the phase space. In general, 
more slowly- varying functions will reveal only the broad fea- 
tures of the phase space, while fast-variation and localized 
features can reveal smaller-scale dynamics. 

As an example, we consider functions from the Haar 
wavelet basis in 2D. Wavelets are a multi-scale functional 
family used in frequency decomposition and multi-resolution 
analysis [35]. In ID, the Haar wavelet basis is constructed 
from the mother wavelet ipoo, defined by: 

r 1 if xe [0,|] 
ip 00 (x) = \ 1 if aj€]|,l[ 

[ if ieR\[0,l] 

A general ID wavelet ipy is constructed from this one by 
the appropriate transformations, summarized as: 

^(x) :=?/ 2 i; 0Q (?x~2- l j), 

and, as it can be shown, the functional family {ipij}i,jez 
is a basis for L 2 (R). The 2D wavelets are constructed as 
products of ID wavelets (in case 2D functional family is 
to be a basis for L 2 (IR 2 ), the procedure also involves cross- 
products of ID wavelets with ID scaling functions; we shall 
however skip these technical details, and refer the reader to 
[28] and references therein). For the purposes of our study, 
we construct and use the 2D wavelet function called W: 

W = 2 ^ £n,m=l <» ^3ra = (g) 

= El, m =i <M8x - f ) ® VW8* - f ), 1 ' 

whose graph is shown in Fig. [2^,. Clearly, W can be con- 
structed from the 2D wavelets basis. Note that W is not 
continuous and thus does not strictly correspond to theory 
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Figure 1: Single-function color-plots of time averages for the standard map Eq.([7| (left), and the corresponding phase space 
portraits (right) done as 100 iterations of 11 x 11 random trajectories picked from [0, l] 2 , for the function / = cos(27ry). 
Top row: e = 0, middle: e = 0.09, bottom: e — 0.18. The grid of 800 x 800 initial points was used, and the dynamics was 
run for ig na i = 30000 iterations. 



in pp. However, this can be readily remedied using com- 
pactly supported continuous wavelets. 

In Figs.[2]D, c & d we show time averages of W under the 
standard map dynamics for various values of e- values. It 
is again easy to recognize the known standard map features 
visible through diverse coloration of the phase space regions. 
As opposed to the Fourier functions case, wavelet functions 
have "sharp edges", which is why these plots have more 
drastic coloration changes among closely located invariant 
sets (compare FigJ^D to the middle plot in Fig. [I] regarding 
e = 0.09, and Fig.[2p to the last plot in Fig. [T] regarding e — 
0.18). This property can be used for "zooming" (see Section 
[5]), e.g. in the case when a specific phase space sub- region 
is to be examined. It is interesting to note that the chaotic 
region for e = 0.16 (cf. Fig.llb) is not uniformly colored 
despite e being above the chaotic transition value - as the 
transport throughout the chaotic region is still very slow 
in the area around the last broken KAM curve, a different 



coloration occurs there. This type of chaotic transport in 
Hamiltonian maps can be approximated by the Markov tree 
model [231301131]. 

In order to investigate the phase space structure at larger 
(smaller) scale, one takes smaller (bigger) Fourier /wavelet 
frequency. In the context of standard map it is convenient 
to use a slightly bigger frequency for y coordinate - since the 
transport in x direction is much faster, the function averages 
out to zero rather quickly in this direction. Time-average 
plots for several Fourier basis functions with increasing fre- 
quency are shown in Fig.[3j Note the relationship between 
the used frequency and the scale of the best-visualized phase 
space details. 

While the Fourier functions give smoothly colored plots, 
wavelets produce sharper and more detailed plots with a 
better distinction between the independent invariant sets. 
However, different invariant sets still happen to be assigned 
the same colors. To remedy this, we use the full power of 
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Figure 2: The graph of the wavelet W as defined by Eq. ^ in (a). Its time average under the dynamics of the standard 
map Eq. Q for e = 0.09 in (b), e = 0.16 in (c) and e = 0.18 in (d). Grid: 800 x 800 and i final = 50000 iterations. 



the ergodic partition concept in Section [5j 

4 The Convergence Properties 

It is well known that time-averages of functions under dy- 
namics of measure-preserving maps can converge arbitrarily 
slowly (see [3T] for discussion and related results). However, 
for specific continuous functions bounds on rate of conver- 
gence are computable [32] , Bounds that depend only on 
the maximum absolute value of a continuous function can 
be obtained [32], but in our context it is relevant to study 
how convergence properties depend on the type of trajec- 
tory of /. Thus, in this Section we study numerically the 
convergence properties of time averages in order to estimate 
the precision of the values obtained and relate them to the 
different types of orbits. For simplicity we use only Fourier 
functions (the results for wavelets are similar) . Consider the 
i-th partial time average for a function / given by: 

t-i 

/*(z ,lto) = ^/(T^o.Ito)), (9) 

fc=0 

with lim^oo /'(xo, Uo) = f*{x$, yo) (which we assume exists 
for all grid-points (xo,j/o))- The difference: 

A(t) = \r(x ,yo)-f(x ,y )\ (10) 



is a sequence whose asymptotic behavior is to be studied 
in relation to the initial point (xa,yo) and the e-value, de- 
pending on the phase space regions with different dynami- 
cal behaviors. We consider the function / = cos(27ry) (cf. 
Fig.[lj), define /* = /* for t = 10 8 and consider the first 10 6 
iterations. 

The Regular Region. For all the points on the regular tra- 
jectories the time averages of continuous functions converge 
with the error decreasing as § , with the constant a given by 
the trajectory properties, as shown in [19] . The result ap- 
plies uniformly to all the regular (periodic or quasi-periodic) 
orbits, regardless of the e-value and the choice of function. 
Given that a ~ 0(1) this allows a rather precise estimation 
of the final precision of the obtained time average value, in 
relation to the total number of iterations computed f final- ^ 
typical convergence plot for the case of a regular trajectory 
is reported in Fig.[4k. 

Strongly Chaotic Region. In the case of strong mixing 
the fluctuations of the time average decrease as ^ [T9] . 
Our findings indicate this rate asymptotically approaches 
t a regime with a > — |; the a = — | result was obtained 
only in the case of very chaotic orbits (large e). A typical 
convergence plot for a chaotic trajectory is shown in Fig.[4]D. 
Given the mixed phase space of map Eq. (TTb, the plot ex- 
hibits a very irregular behavior, which is bounded by the 
t a regime away from the transients. Hence, the error in the 
case of chaotic time averages can be estimated from below 
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Figure 3: Color-plots for various time averaged functions for the standard map Eq.Q with e = 0.12. More detail of the 
dynamics (higher order resonant islands) are revealed when the frequency of the function is increased. Grid: 800 x 800 
and fgrini = 30000 iterations. 



by V Vernal anc ^ improved by the increase of the e- value. 

Weakly Chaotic Region. Below the chaotic transition for 
Eq- ( £ ~ °- 15 )> the chaotic regions are localized in the 
phase space around the hyperbolic fixed points, and char- 
acterized by a weak chaos with trajectories slowly diffusing 
through the region |33j . The convergence of the time av- 
erages in this region is extremely slow and irregular, and 
cannot be in general bounded by any asymptotic slope of t a 
type. A typical convergence plot for this region is reported 
in Fig.|4j:. Note that despite plot showing less irregular oscil- 
lations than in the strongly chaotic case, it barely decreases 
and cannot be fitted with a determined slope. This conver- 
gence pattern is consistently present in all the weak chaos 
trajectories, and it improves only with the increase of e- 
values. It is therefore hard to estimate the final error in this 
case, unless the dynamics is run for excessively long times. 

Our visualization method appears to suit better the cases 
of either regular or strongly chaotic behavior where the 
"speed" of filling the invariant set is relatively high. How- 



ever, we find that even in the case when weak chaotic be- 
havior is present, we obtain a good representation of the 
structure of the phase space, given that regular and strongly 
chaotic orbits fill a large portion of the phase space. In view 
of this, one can set the total number of iterations according 
to the precision rates of the strongly chaotic zone. We thus 
typically take for standard map % na j ~ O(10 4 ) iterations: 
this sets the precision of strongly chaotic case to O(10 -2 ) 
(with even better precision for the regular case), which is 
enough given that the considered functions have values in 
the [—1,1] interval. Note also that the convergence prop- 
erties show a certain pattern in relation to the trajectory 
type: it would be possible to characterize the nature of a 
trajectory by looking at this patterns for different functions 
(similarly as done in [H]). In Figs.[9|fc[l5| we examine the 
convergence slopes for more functions, and averaged over 
and ensemble of trajectories. 
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(a) (b) (c) 

Figure 4: The convergence plots of the time averages for the standard map dynamics. The function / = cos(27ry) is 
considered, and /* is taken to be f 10 . (a): a regular orbit for e = 0.09 (cf. Fig.[l]) for the initial point (0.5,0.4) fitted 
with the slope of -1; (b): a strongly chaotic orbit for e = 0.18 (cf. Fig. [I]) for the initial point (0.02,0.02) fitted with the 
slope of — |; (c) a weakly chaotic orbit for s = 0.09 (cf. Fig.nj) for the initial point (0.01, 0.01). 



5 Ergodic Quotient Space and Clus- 
tering Methods 

In this Section we present the concept of the ergodic quotient 
space and an algorithm based on time averages of functions 
that helps us understand its topology. 

A discrete-time measure-preserving map x' = Tx on 
a compact phase space A decomposes the subset of the 
phase space £ on which time averages of all continuous func- 
tions exist into ergodic sets [TJ The ergodic quotient 
space Q e is the space where each ergodic set is mapped 
into a single point, and, additionally, the complement of 
£ is mapped into a point. Consider the space S of all 
the infinite sequences indexed by non-negative integers S = 
{do, ai, a 2 , ...}. Consider also a basis i 6 N for L? on A. 
Denote the image of q : £ -> S, g(x) = {xs<= (x) , /J (x) , /* (x) , . 
by Qe, where xs c is the indicator function on set £ c and 
call the map q the ergodic quotient map. For convenience 
we extend the definition of /* from £ to the whole set A by 
setting /* = on £ c . Then, clearly, 

Q e = Q* e u {1,0,0,...}. 

Now note that finite-dimensional Euclidean space embed- 
dings of projections of Q e to a space of finite sequences la- 
beled by ii , ijv can be obtained by iV-tuples {/£,..., /* } (x 
We call such embeddings Mesochronic Scatter Plots and 
study them numerically in this Section. But before that, 
let us give an example for ergodic quotient space of the dis- 
crete dynamical system presented in Eq. Q with e = 0. 
In that case, fixing an irrational y, the time average of any 
function is constant, for any initial x. Thus, the whole circle 
y = const, is mapped into a single point. For rational values 
of y = | ( | is an irreducible fraction), where the invariant 
circle at that y is filled with periodic orbits, the length of 
the interval parameterizing distinct periodic orbits is K It 
is then easy to see that Q e is a "rational comb", consisting 
of a straight line with an interval of length | attached at 
every rational y = |. 



5.1 2-dimensional MSP embedding 

To provide a graphical representation of Q e , we consider 
Fourier basis functions and begin by the case of N = 2 
basis functions. Time averages /* and /| are computed for 
every grid-point defining the correspondence between the 
grid-points (xo,j/o) & n d the time average vector {/*, /I}: 

(x ,yo) € grid — ► (fi(x a ,yo)J^(x ,y )) € [-1,1] 2 . 

The MSP is obtained by plotting all the vectors (/* , /|) on 
the square [— 1, l] 2 . We use the grid of 300 x 300 points and 
run the dynamics for 30000 iterations. 

In Fig.[5]we show a sequence of 2D MSP for the functions 
fi = sin(27ry) and f% — cos(127ra;) cos(27r?/) with increasing 
e- value. In the first plot showing the case of e = 0, for irra- 
tional y, the time average of f\ is sin(27ry), while the time 

■aVerage of fs is 0. This explains the horizontal line in the 
e = plot. For rational y the time average of /3 is only non- 
zero provided y — 0, |, |, |, | and y = | (in other words, 
the function "resonates" with the map dynamics and pro- 
duces non-zero values of time averages only for those y's). 
These values contribute the vertical lines in the e = plot. 
However, some of these vertical lines overlap in this projec- 
tion, which will lead us to consider three-dimensional pro- 
jections later. The side lines will become more evident with 

I increasing e, due to the appearance of resonance zones and 
the associated new families of quasiperiodic orbits. With 
further increase of e, the side lines disappear and a central 
scattered region appears, leading to the chaotic transition 
at the familiar value close to e « 0.154. 

The idea is further illustrated in Fig. [6] where a MSP in- 
volving two functions is investigated with reference to the 
corresponding time averages plots and the phase space por- 
trait. As already stated, each single point in the MSP has 
the ^-coordinate equal to fi(x ,y ) and the y-coordinate 
equal to (xq, yo) for some grid-point (xo, yo)- As indicated 
in the figure, long branches (curves) represent the families 
of periodic islands around elliptic fixed points, while the 
irregular clouds amount for localized chaos around the hy- 
perbolic fixed points. Secondary chaotic zones appearing 
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Figure 5: Sequence of 2D MSPs for standard map Eq. ([7]), with functions /i = sin(27ry), /3 = cos(127ra;) cos(27ry). The 
value of e is indicated in each plot. Each time average was obtained on a 300 x 300 grid, for 30000 iterations. 



around second-order hyperbolic points are also visible, to- 
gether with the secondary families of periodic orbits. 

To investigate this further we show the same MSP for 
the same time averages, zoomed to the phase space region 
[0.6, 0.9] x [0.6, 0.9] in Fig.[7j This MSP can readily be rec- 
ognized as a part of the MSP from Fig. [6] Again, we see the 
interplay between long curved lines and irregular clouds, 
representing regular and chaotic regions respectively. Note 
that many more secondary periodic families are visible in 
this plot due to the improved resolution (zoom), captur- 
ing the scale-invariant fractal nature of the standard map's 
phase space. 

5.2 3-dimensional MSP embedding 

The two-dimensional projections presented in the previous 
Section have the unpleasant feature of self-intersection. Three- 
dimensional embeddings resolve this issue for two dimen- 
sional maps. We will see later that higher dimensional em- 
beddings are needed to avoid intersections in higher dimen- 
sional maps. 

We set N — 3 and consider the MSPs done with three 
linearly independent functions. Using the same grid and to- 
tal iterations as previously, we consider the correspondence: 

(x ,y Q ) — > (ft(x ,yo), /lOcs/o), fSi x o,yo)) e [-M] 3 , 



with f 2 = cos(27ry). In Fig.[8] we show nine MSPs obtained 
for the three functions and increasing e-values. We mon- 
itor the phase space structure evolution as e is changed, 
in terms of geometric complexity evolution of the MSP's 
structure, which is much easier to do within this 3D embed- 
ding. The size of resonance zone emanating from circle of 
period-6 orbits is substantially shrunk compared with that 
for period-3, period-2 and period-1, even for small, e = 0.01 
perturbation (middle figure of top row in Fig.[8|. Note the 
changes in branches and development of higher order peri- 
odic islands, with increased e with localized chaotic zones 
around hyperbolic periodic orbits. These secondary islands 
have resonance zones of their own that merge and enable the 
chaotic transition, where a single large chaotic zone enables 
trajectories to pass around the torus in both directions. The 
chaotic regions are visible as thickened scatter in the plots. 
The chaotic transition occurring for e ~ 0.154 can here be 
seen as a merging of localized chaotic zones that propagate 
along the branches into a single connected chaotic zone, vis- 
ible as a single cluster for the value e = 0.18 in Fig. [8] Also 
visible in the e = 0.18 plot are the side lines correspond- 
ing to the remaining islands around y = 0, |, | (see bottom 
plot in Fig. [I]). As expected, in the strongly chaotic regime 
all the time average vectors are localized in a single cluster 
that shrinks in size with further increase of e. In the case 
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Figure 6: Two-function MSP (up left) for the standard map Eq. ^ for e = 0.12 (cf. Fig. [5]). Time-averaged functions 
are /i = sin(27ry) (down left) and fs = cos(127ra) cos(27ry) (up right). Phase space dynamical regions are indicated, in 
correspondence with the MSP parts and the dynamical regions in phase space portrait (down right). 




Figure 7: A zoomed part of the MSP (up left) from Fig. [6] for the phase space region [0.6,0.9] x [0.6,0.9] (for s = 0.12, 
fx = sin(27ry) (down left) and fa = cos(127ra;) cos(27r?/) (up right)). Dynamical regions and their corresponding MSP parts 
are indicated are related to the phase space portrait (down right). 



of ergodic behavior one expects all the time average vectors 
to shrink to a single point in [— 1, l] 3 space. 

Adding a function to the MSP, and thus increasing the 
embedding dimension clearly improves the representation 
of the phase space structure. Furthermore, the Fig.|H]shows 
the sufficient embedding dimension for the standard map's 
MSPs to be three: it is in three dimensions where the fami- 
lies of regular orbits can be fully represented without inter- 



section. 

Finally, in Fig. [9] we examine the time-evolution of the 
time average vectors for e = 0.18, shown in the last plot in 
Fig.|] Time averages of the same functions are computed 
for various final iteration- values ifj na i, and for each grid 
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Figure 8: Sequence of 3D MSPs for the standard map Eq. Q, with functions /i = sin(27ry), /2 = cos(27ry) and fa — 
cos(127rcc) cos(27ry). The value of e is indicated in each plot. Each time average was obtained on a 300 x 300 grid, for 
30000 iterations. 



point (xo,yo) the norm of time average vector 

|f*(*o,yo)| = v/(/!(^o,yo)) 2 + (/J(zo,j/o)) 2 + (fl(x ,y )) 2 

( n ) 

is considered. The distribution of values of |f*| is shown in 
Fig. [9^i as function of time t. While for large |f'|-values we 
see a quick convergence into a final profile, for small |f*|- 
values the distribution-profile seems to be slowly evolving 
towards a sharp single-peaked (delta) distribution. Clearly, 
the former corresponds to the regular orbits that are faster 
to converge to final /*-values, while the latter corresponds 
to the slowly converging chaotic orbits (cf. Section [4J. We 
illustrate this further be generalizing the Eq.[H)]for the case 
of more functions into: 



A(t) = ||f*(x , W ))|-|f t (a:o,W))|| 



(12) 



Thus, A(£) measures how close is the norm of the partial 
time average vector |f'| to its limit value |f*|. In Fig.^ we 
show the distribution of A(t) for all grid-points in function of 
time t, obtained by taking |f*| = |f *=i-28x 10 |_ rp w0 g r0U p S 
of peaks that travel to zero (lnA(<) — > —00) and statisti- 
cally maintain their shapes can be readily recognized. The 
structured group consisting of few smaller peaks has much 



smaller In A(t)-values than the single structureless peak lo- 
cated around In A(t) ~ —3. Also, the former group of peaks 
seems to travel about twice faster towards zero than the lat- 
ter single peak. This again corresponds to the difference be- 
tween regular orbits (structured peak group) and the chaotic 
orbits (the single structureless peak) - the former converge 
to zero with a rate of i -1 uniformly for every point, while 
the latter on average converge to zero with the rate of t~ 5 . 

5.3 Visualization of the ergodic partition 
via clustering 

Following the investigation of the MSPs we construct a sim- 
ple algorithm for approximation of the ergodic partition 
and its graphical phase space visualization. Consider an 
N- function MSP contained in [—1, ll 



liV. 



step 1 divide the TV-cube [—1,1] into L cells dividing 



each axis into L segments as illustrated in Fig. 10 1 
(for N = 2 and L = 10), and consider the distribution 
of time average vectors around the cells 

step 2 disregard the cells that contain no time average vec- 
tors 
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Figure 9: Time averages for the functions fi — sin(27ry), f% 
various iteration-values for the standard map Eq. ^ on the j 
the time-evolution of the norm of time average vectors |f*| = 
|f'| as function of time; (b): distributions of A(t) = ||f*| — |f 



= cos(2tt?/) and fs = cos(127ra;) cos(27ry) are computed for 
»ri d of 400 x 400 with e = 0.18 (cf. last plot in Fig.[8]), and 
V Hi) 2 + (il) 2 + {fi) 2 1S considered, (a): distributions of 
: || as function of time, done by taking |f*| = |f*=i-28xio |_ 



step 3 assign a color to every remaining cell, therefore as- 
signing a color to every time average vector 

step 4 observe the grid-points corresponding to time av- 
erage vectors sharing a cell/color: they define an iV- 
order approximation of an ergodic set 

step 5 color the phase space by coloring each grid-point 
with the color assigned to it 

Note that this is a generalization of the single-function col- 
oring scheme with a difference that now the scheme is regu- 
lated by adjusting the cell division and optimizing it accord- 
ing to the structure of the MSP. The multi-dimensionality 
of the color-assigning rule allows for higher differentiation 
among the invariant sets. 

For simplicity we start again with the case of two func- 
tions: one-dimensional lines intersect only at points and this 
feature does not appear to perturb the phase space represen- 
tation much. We examine the two-function MSP showing it 
in Fig.[T0^, with 10 x 10 cells division. The corresponding 
approximation of the ergodic partition with colors between 
blue and red randomly and uniformly assigned to non-empty 
cells is shown in Fig.|l0p. A better overall clarity of the in- 
variant set structure is obtained both at the global and local 
level (within the approximation precision, which is also in- 
fluenced by a limited number of available colors). Note that 
the visibility can be enhanced by coloring nearby invariant 
sets with different colors which is attained by randomizing 
the color-assignment for the non-empty cells. The number 



of visualized sets increases with increase of L (Figs. 10 3 & c) 



but the color differentiation gets poorer, as the number of 
available (visible) colors remains limited (not only by the 
software, but also by the human eye recognition). The op- 
timal value of L (regulating the number of cells) is to be 
set according to the visualization requirements, taking into 
account the relationship between color differentiation vs. 
number of visible invariant sets. For the two-function case 



examined in Fig. |10[ it appears the optimal L is around 50 
(corresponding to Fig.[To};). Too small L is underusing the 
MSP as it colors too many different invariant sets uniformly 
(Fig. |To]b) , while for too large values of L the lack of colors 
brings the same problem, in addition to a very non-uniform 
coloration of the chaotic region (Fig.[To}l). 

In Fig.[TT] on the left we show an example of ergodic 
partition approximation constructed from a three-function 



MSP, using the functions from Fig. 10 case as the first two. 
Note that for L — 50 we obtain a better quality than previ- 
ously. Finally, in Fig.[TT]on the right we add another func- 
tion and show a four-function approximation obtained for 
the optimal L value of L = 50. 

As noted earlier, we have chosen functions involving dif- 
ferent frequencies, thus visualizing global and local phase 
space features simultaneously. Fig.[TT]revcals high-resolution 
approximations to the ergodic partitions for e = 0.12, pro- 
ducing good approximations to the invariant set structure of 
the standard map's phase space. Note that both pictures in- 
deed visualize details at all scales, with the plot on the right 
being somewhat sharper. In construction of these plots we 
sought to improve the cell division by having a relatively 
uniform number of time average vectors within each cell, in 
relation to the available colors. Due to a particular choice of 



functions, in the right plot on Fig. 11 we managed to obtain 
good coloration for L = 50 (very large L-value considering 
the number of functions involved) , creating the optimal ap- 
proximation for the four-function case. Given the number 
of invariant sets visualized, different realizations of random 
colors assignments make very little difference in the overall 
picture. 

The limited number of available colors still makes some 
different ergodic sets appear in the same color; this problem 
can be partially overcome by optimizing the color-assigning 
rule. Instead of assigning a random color to each non-empty 
cell one could create an assigning algorithm that would be 
optimized in relation to the particular case studied. An- 
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Figure 10: Two-function approximation of the ergodic partition, (a): two-function MSP for the standard map Eq. |7| for 
e = 0.12, done using cos(27ry) and cos(27rx) cos(27ry), for a grid 800 x 800 and with % na j = 30000. Three approximations 
of the ergodic partition done for different values of cell division L (and constructed using this MSP for clustering), are 
shown in (b) for L = 10, in (c) for L = 50, and in (d) for L = 140. 
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Figure 11: Approximations of the ergodic partition for the standard map Eq. ^ with s = 0.12. Left: three-function 
approximation using two function from Fig. 10 and the function sin(47rx) sin(47ry). Right: four-function approximation 



using these three functions in addition to sin(107ra) sin(107ry). The grid 800 x 800 is used for all time averages and the 
cell division for both approximations is L = 50. 



other problem arises in relation to the ergodic zone in the 
phase space: given that its time average vectors are more 
diffused then the regular orbits' ones (cf. Fig. [8]), they set 
the lower bound to the size of the cells. This is why the 
ergodic zone in the Fig. 10 1 appears non-uniformly colored. 



In the context of the standard map, it is convenient to pick a 
smaller L for low e-values in order to obtain a better focus 
on the nested invariant curves, while for larger er-values a 
bigger L allows to include the whole chaotic zone in a single 
cell/color. Also, a bigger number of functions allows more 
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flexibility for the L-value as the underlying MSP differen- 
tiates better among the invariant sets. Moreover, a better 
cell division scheme would not employ the simple cubical 
division described here, but a more sophisticated algorithm 
optimizing (for instance) the number of time average vec- 
tors per cell with respect to their distance in the time aver- 
age space. One could also seek to adjust the shape of cell 
according to the properties of the MSP instead of just us- 
ing the simple cubical ones. Furthermore, the problem of 
limited number of colors could be tackled by employing a 
specific graphically-oriented visualization software allowing 
more flexibility in terms of choosing or adjusting colors and 
their tones. Also, a cell's color could be determined in re- 
lation to the number of time average vectors contained in 
it, or contained in the neighboring cells, thus differentiating 
better among various invariant sets. 

Finally, in the case of measure-preserving maps with a 
given maximum MSP embedding dimension (three in this 
case, cf. Fig. [8]), one could seek to parameterize the MSP 
obtaining a continuous coloration scheme that would include 
all the invariant sets. Still, this procedure would yield a non- 
uniform number of time average vectors per cell, allowing 
for further improvements. For the purpose of this study 
however, we limit ourselves to the simple and illustrative 
algorithm just exposed. 

6 The Froeschle Map 

As our first higher dimensional example we consider the 4D, 
measure-preserving Froeschle map [22] that consists of two 
standard maps with a symplectic coupling: 

x\ = x\ + yi + E\ sin(27ra:i) + r\ sin(27ra;i + Ikx-i) 
y[ = j/i + £i sin(27rxi) + r7sin(27ra;i + 27TX2) 



y'2 



X2 + 2/2 + £2 sin(27rx 2 ) + ?7sin(27rxi + 27ra 2 ) 
2/2 + £2 sin(27rx2) + rjS\n.{2iTXi + 2-kx-i) 



(13) 

where (xi, yi, x 2 , J/2) € [0, l] 4 - We set E\ = e 2 = £ and re- 
duce our investigation to the case of two identical interacting 
standard maps. This map can be related to the standard 
map in two ways: for 77 = we have a system of two uncou- 
pled standard maps, while for e = one can introduce the 
new variables: 

U\ = X\ + X2 

v\=yi + 2/2 

U2 = Xi — X 2 

v 2 =yi- 2/2 

which reduce the system to two maps: a standard map in co- 
ordinates (ui,vi) with the parameter 2rj, and an integrablc 
twist map in [u2,v 2 ). Given this structure, an interesting 



way of varying parameters of the map Eq. ( 13 ) is by letting 
e = 2-q. We slice the 4D phase space by a 2D section with 
fixed (£2,2/2) = (0,0), and consider a grid of 500 x 500 ini- 
tial grid-points in (x\, 2/i)-space, plotting the time averages 
of a single function f 2 = 2 cos(27rj/i) + cos(27ryi) cos(27r?/2) + 
cos(127ra2) + cos(127ra;i) obtained by running the dynamics 
for 200000 iterations. The results are shown in Fig.[l2j Note 
the departure from the known standard map's phase space 



due to increasing interaction between the maps. For small s 
values, the phase space maintains its regular structure with 
only small and localized chaos similarly to standard map for 
small s. The coloring here indicates two dimensional inter- 
sections of the real 4D invariant sets with the selected 2D 
phase space section. It appears that the global chaotic tran- 
sition in the sense of vertical transport through the phase 
space occurs around e = 2-q = 0.05. 

The structure of Q e for the Froeschle map for e = r\ = 
is clearly a product of two "rational combs" discussed 
earlier. This product is topologically a torus to which: 

• a line of length ^ is attached at every point (2/1,2/2) 
such that one of the j/j's is irrational and the other 
rational, | (with ^ an irreducible fraction), 

• a rectangle of sides qi , q 2 is attached at every point 



(2/1 = ^,2/2 
fractions. 



— ), where — and — are irreducible 



Since two rectangles generically do not intersect in a 5 di- 
mensional Euclidean space, the appropriate embedding di- 
mension for the Froeschle map is five. To avoid self-intersections 
in the embedding, we would thus need to consider at least 
five independent functions. This information is important 
when choosing the number of functions that are used in clus- 
tering and approximation of the ergodic partition. However, 
since we can not graphically represent the 5D space results, 
we consider here three-function MSP visualized in 3D em- 
bedding. An example with e = r\ = is shown in Fig. 13 
obtained for a 4D grid 13 x 13 x 13 x 13 initial points as 
follows: 

• the lattice-grid is set in full 4D phase space, and iter- 
ated using the Froeschle map Eq. ( 13 ) with e = r] = 1 
for 10 iterations, in order to randomize the points 
(starting from the uniform lattice-grid in the case of 
e = 77 = might create trajectories that strongly over- 
lap 

• from thus obtained 13x13x 13x13 = 28561 initial 4D 
points, set 13 x 13 x 13 = 2197 to have (j/i =0,y 2 = 0), 
another 13 x 13 x 13 to have (2/1 = \,y% = \) and 
so on, depending on how many resonances are to be 
visualized 

• run the dynamics for tg na } = 100000 iteration starting 
from thus created ensemble of initial 4D points 

• as the result, majority of points will account for quasiperi- 
odic orbits, while the selected points will visualize the 
chosen resonances 



Note that in the Fig. 13 all the mentioned phase space fea- 
tures are visualized as expected. The chosen resonances for 
this case are (2/1 = 0,y 2 = 0), (2/1 = 5,2/2 = §), (2/1 = 
|,2/2 = |), and (2/1 = §,2/2 = f )■ 

We examine the three-function MSPs for the Froeschle 
map with e — 2 77 > 0, by employing 4D grid of 12 x 12 x 
12 x 12 initial grid-points and the same number of total 
iterations. The ensemble of initial points with selected reso- 
nance points is constructed as above. We consider structural 
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Figure 12: Time averages of a single function f 2 = 2cos(27ryi) + cos(27ryi) cos(27n/2) + cos(127ra;2) + cos(127rai) under 



the dynamics of the Froeschle map Eq. (13), computed for a 500 x 500 grid of initial conditions taken on the phase space 



section (x 2 = 0,y2 = 0). The dynamics was run for 2 x 10 5 iterations. The values of e and i], linked by e = 2rj, are 
indicated in each plot. 



changes in the MSPs by changing the e-value, as shown in 
Fig. 14 As the coupling intensity grows, the structure from 



the previous figure is destroyed, in a way similar to what 
observed for the standard map (cf. Fig.[8|. In particular, 
note the different mechanisms of destruction of resonances 
and irrational orbits. The MSP's structure reports a given 
level of regularity in the phase space persisting for a certain 
range of coupling parameter strengths. This is a manifes- 
tation of K AM/resonance zone nature of this map, which 
(similarly to the standard map) maintains some invariant 
tori until the parameters exceed certain thresholds. From 
the MSP, we conclude that global merging of resonances oc- 
curs for value of e = 2-q between 0.05 and 0.06 (cf. phase 



space cross sections visualized in Fig. 12 | 



We conclude the Section by observing that MSP analysis 
as employed here allows investigations of systems much more 
complex than the standard map. As pointed out previously, 
with an appropriate choice of functions one can reduce the 



complexity of a given system to geometrical features of the 
MSP, capturing the key dynamical details of the system in 
the form of MSP structure. Moreover, various changes in 
system's properties can be monitored this way (cf. Fig. 14) 
by observing the geometrical evolution of the MSP. A fur- 
ther application of this technique might be in the study of 
invariant sets of even higher dimensional dynamical systems, 
like the coupled maps on networks [25j [26] and recently dis- 
covered maps that mimic quantum chaos [35 . 



7 Extended standard map 

As a second higher dimensional example we consider the 3D 
extended standard map [3], that represents a generalization 
of the classical standard map. It is a volume-preserving 



15 



e=0.01, 11 = 0.005 



e=0.02, r| = 0.01 



e=0.03, i]=0.015 






£ = 0.04, r| = 0.02 



e=0.05, ii = 0.025 



e = 0.06, ri = 0.03 



„_m - 






£ = 0.07, n = 0.035 



e=0.08, \\-0.M 



£=0.09, ii = 0.045 



mm 




Figure 14: Sequence of three-function MSP for the Froeschle map Eq. (13). The grid of 12 x 12 x 12 x 12 initial points 



is constructed as for the Fig. [13] The functions used are: /i = 2sin(27ryi) + sin(27ryi) cos(27rj/ 2 ) + cos(127rxi), f% = 
2 cos(27ryi) + cos(27n/i) cos(2-7r?/2) + cos(127ra;2) -I- cos(127rxi) and /a = sin(2-7ry2) + cos(127r:C2), and the dynamics was run 
for 10 5 iterations. The respective e = 2rj values are indicated in each plot. 



action-action-angle map defined as: 

x' = x + e sin(27rz) + S sin(27r?/) 

y' — y + e sin(27rz) 

z' = z + x + e sin(27rz) + 6 sm(2i:y) 



[mod 1] 

[mod 1] (14) 
[mod 1] 



with the values in [0, l] 3 . Its physical origin and the ana- 
lytical properties are investigated in [3] . For 6 = the map 
takes the form: 



x = x + e sin(27rz) [modi] 
y 1 = y + e sin(27rz) [mod 1] 

z' = z + x + e sin(27rz) [mod 1] 



(15) 



which keeps the planes y — x = const, invariant under the 
dynamics, and reduces to the standard map in x and z co- 
ordinates (while y behaves like another action coordinate). 
This family of standard maps on diagonal planes is however 
broken for 5 > as the transport is allowed between the 
diagonal planes. 

It was conjectured in [3j that the extended standard 



map Eq. ( 14 ) is ergodic for small positive values of per- 
turbations e and 6. We bring additional evidence to this 



vectors shrink to zero with time-evolution of this map for 
e = 0.01 and 5 = 0.001 (equivalently to what was done 
in Fig.[9|. Time average vectors f* = (Z*,/!)/!) ar< 3 corn- 



claim in Fig. 15 where we examine the way time average 



puted for extended standard map Eq. ( 14 1 with e = 0.01 
and 6 = 0.001. The grid of 50 x 50 x 50 initial points 
was used, with randomized points as in the previous Sec- 
tion. The functions /i = sin(27nc),/2 = cos(27ry) and / 3 = 
sin(47ry) cos(47ra;) cos(127rz) are considered. In Fig. 15 1 we 
show the distributions of time average vector norms |P| de- 
fined as in Eq.JTl] in function of time t (cf. Fig. [9^). It 
is interesting that at t = 400000, the structure of reso- 
nances is present due to a quasi-two-dimensional nature of 
the map causes a multi-modal distribution of time averages. 
This multi-modality disappears with the higher number of 
iterates and the distribution assumes an exponential shape, 
with variance and mean tending to zero. In Fig. |15b the 
mean values of distributions are shown in function of number 
of iterations t, and fitted with the slope of -0.21. This sug- 
gests an average convergence rate of t -0 ' 21 , slower than for 
the case of strong chaos characterized by t~°- 5 (cf. Fig.]^). 

This seems to confirm the mentioned ergodic hypothesis 
from the computational prospective. It also agrees with the 
result proved in [3] stating that no invariant two-dimensional 
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Figure 15: Time average vectors f* = (/f, /Ji /J) are computed for extended standard map Eq. (14 1 for the functions 
fx = sin(27ra;), f% = cos(27r?/) and fs = sin(47ry) cos(47r:r) cos(127tz) with e = 0.01 and 6 = 0.001 on the grid of 50 x 50 x 50 
random initial points for various numbers of iterations, (a): time-evolution of the distribution of time average vector 
norms |f*|; (b): the distribution's mean value as function of time t, fitted with the slope of -0.21. 



8 Conclusions 




We presented the computational realization and theoretical 
extension of an invariant set visualization method based on 
ergodic partition theory suggested in [2] . We defined the er- 
godic quotient space obtained by associating an ergodic set 
with a point. Embeddings of the ergodic quotient space into 
Euclidean space, called Mesochronic Scatter Plots (MSP's) 
were realized using time averages of observables on the phase 
space. The time averages were computed and visualized us- 
ing a coloring scheme that we named a Mesochronic Plot 
for a variety of measure-preserving maps. The time average 
convergence issues have been considered. A simple algo- 
rithm for approximation of the ergodic partition was de- 
veloped from multi-functional MSP's by dividing the time 
average vectors space into cubical cells. Approximation of 
ergodic partition structure in the phase space was shown 
for various numbers of functions. By studying the standard 
map with known properties, we were able to confirm the vi- 
sualization results within the limits of numerical precision. 
The extent of the method's applicability was illustrated on 
4D Froeschle map and 3D extended standard map, giving 
new insights into dynamical structure of these systems. Er- 
godic invariant sets in higher dimensional systems can be 
visualized using our method by obtaining their intersections 
with two-dimensional surfaces of choice. 

In the paper to follow [4] we show how periodic sets 
surface persists in this map for any positive perturbation and resonances can be graphically visualized according to 
value, implying that the map allows a global transport throughiheir periodicity and the phase space structure, by the use 



Figure 13 



Eq.(13| 
13 



Three-function MSP for the 
with e = r\ = 0. The grid 
13 initial points is constructed to 



Froeschle 
of 13 x 
visualize 
The dynam- 

run for 10 5 iterations. The used functions are: 



map 
13 x 
both 



resonances and irregular orbits (see text). 



ICS if 

/i = 2sin(27ryi) + sin(2-7ryi) cos(27n/ 2 ) + cos(127rxi), f 2 = 
2 cos(27ryi) + cos(27ryi) cos(27rj/2) + cos(127ra;2) +cos(127rxi 
and /3 = sin(2-7ry2) + cos(127ra;2). 



out the phase space at a small non-zero perturbation. This 
consideration demonstrates our method to be useful in the 
context of numerical investigations related to ergodic prop- 
erties of dynamical systems as the one discussed above. 



of harmonic time averages that extend the concept of time 
average described here. Both the method presented here and 
the method in the follow-up paper are related to eigenspace 
structure of the Koopman operator. 

The further improvement of the method can be obtained 
by optimization of clustering techniques beyond the simple 
cell division exposed here 36 . The selection of optimal 
functions for embedding is a wide-open question. Moreover, 
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a more detailed geometric analysis of ergodic quotient space 
might yield additional insights into the dynamics, and indi- 
cate a way to construct better algorithms for clustering of 
time average vectors and ergodic partition approximation. 
The integrable maps studied here have interesting mathe- 
matical structure that is non-smooth but still in some sense 
regular. For example, the unperturbed standard map has 
ergodic quotient space structure of a unit circle with an in- 
terval of length - attached at every rational point | . 

The method can also be applied to the continuous-time 
systems for which the ergodic theory results are equally 
valid, and hence the results shown here apply directly. Of 
course, the computation of time averages for a continuous- 
time system is far more numerically demanding. 
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